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We present a theoretical analysis ol dilute gas Bose-Einstein condensates with dipolar atomic 
interactions under rotation in elliptical traps. Working in the Thomas-Fermi limit, we employ the 
classical hydrodynamic equations to first derive the rotating condensate solutions and then consider 
their response to perturbations. We thereby map out the regimes of stability and instability for 
rotating dipolar Bose-Einstein condensates and in the latter case, discuss the possibility of vortex 
lattice formation. We employ our results to propose several novel routes to induce vortex lattice 
formation in a dipolar condensate. 
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I. INTRODUCTION 

The successful Bose-Einstein condensation of 52 Cr 
atoms [TJ |5J [5] realizes for the first time Bose-Einstein 
condensates (BECs) with significant dipole-dipole inter- 
actions. These long-range and anisotropic interactions 
introduce rich physical effects, as well as new opportu- 
nities to control BECs. A basic example is how dipole- 
dipole interactions modify the shape of a trapped BEC. 
In a prolate (elongated) dipolar gas with the dipoles po- 
larised along the long axis the net dipolar interaction is 
attractive, whereas for an oblate (flattened) configura- 
tion with the dipoles aligned along the short axis the net 
dipolar interaction is repulsive. As a result, in compari- 
son to s-wave BECs (which we define as systems in which 
atom-atom scattering is dominated by the s-wave chan- 
nel), a dipolar BEC elongates along the direction of an 
applied polarizing field [H[S]. 

A full theoretical treatment of a trapped BEC involves 
solving the Gross-Pitaevskii equation (GPE) for the con- 
densate wave function (HJ [7] . The non-local nature of the 
mean-field potential describing dipole-dipole interactions 
means that this task is significantly harder for dipolar 
BECs than for s-wave ones. However, in the limit where 
the BEC contains a large number of atoms the problem 
of finding the ground state density profile and low-energy 
dynamics simplifies. In a harmonic trap with oscillator 
length aho = y/h/(mu)), a BEC containing N atoms of 
mass m which have repulsive s-wave interactions charac- 
terized by scattering length a enters the Thomas-Fermi 
(TF) regime for large values of the parameter TVa/ciho 
[B][7]- In the TF regime the zero-point kinetic energy can 
be ignored in comparison to the interaction and trapping 
energies and the Gross-Pitaevskii equation reduces to the 
equations of superfluid hydrodynamics at T = [5J[71[5]. 
When applied to a trapped s-wave BEC these equations 
are known to admit a large class of exact analytic so- 
lutions [9j. The TF approximation can also be applied 
to dipolar BECs [TU]. Although the resulting superfluid 



hydrodynamic equations for a dipolar BEC contain the 
non-local dipolar potential, exact solutions can still be 
found [111 I12j and we make extensive use of them here. 
The calculations in this paper are all made within the 
TF regime. 

Condensates are quantum fluids described by a macro- 
scopic wave function vb = yfpexp[iS], where p is the con- 
densate density and S is the condensate phase. This 
constrains the velocity field v — {Ti/m)V_S to be curl-free 
V x v = 0. In an experiment rotation of the conden- 
sate can be accomplished by applying a rotating ellip- 
tical deformation to the trapping potential [T31 Q3] . At 
low rotation frequencies the elliptical deformation excites 
low-lying collective modes (quadrupolc etc.) with quan- 
tized angular momentum which may be viewed as sur- 
face waves (and which obey V x v = 0). Above a certain 
critical rotation frequency vortices are seen to enter the 
condensate and these satisfy the W_x v = condition by 
having quantized circulation. The hydrodynamic equa- 
tions for a BEC provide a simple and accurate descrip- 
tion of the low-lying collective modes. Furthermore, they 
predict these modes become unstable for certain ranges 
of rotation frequency |151 116) . Comparison with experi- 
ments [I3J E] and full numerical simulations of the GPE 
[T7l fT8j [19] have clearly shown that the instabilities are 
the first step in the entry of vortices into the conden- 
sate and the formation of a vortex lattice. Crucially, 
the hydrodynamic equations give a clear explanation of 
why vortex lattice formation in s-wave BECs was only 
observed to occur at a much greater rotation frequency 
than that at which they become energetically favorable. 
It is only at these higher frequencies that the vortex-free 
condensate becomes dynamically unstable. 

Individual vortices [201 EH [22] and vortex lattices 
[251 IM1 EZ5] in dipolar condensates have already been 
studied theoretically. However, a key question that re- 
mains is how to make such states in the first place. In 
this paper we extend the TF approximation for rotat- 
ing trapped condensates to include dipolar interactions, 
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building on our previous work (5^1127]. Specifically, start- 
ing from the hydrodynamic equations of motion we ob- 
tain the stationary solutions for a condensate in a rotat- 
ing elliptical trap and find when they become dynami- 
cally unstable to perturbations. This enables us to pre- 
dict the regimes of stable and unstable motion of a ro- 
tating dipolar condensate. For a non-dipolar BEC (in 
the TF limit) the transition between stable and unstable 
motion is independent of the interaction strength, and 
depends only on the rotation frequency and trap ellip- 
ticity in the plane perpendicular to the rotation vector 
[T51 [TrJj . We show that for a dipolar BEC it is addition- 
ally dependent on the strength of the dipolar interactions 
and also the axial trapping strength. All of these quan- 
tities are experimentally tunable and this extends the 
routes that can be employed to induce instability. Mean- 
while, the critical rotation frequency at which vortices 
become energetically favorable f2„ is also sensitive to the 
trap geometry and dipolar interactions |21j . and means 
that the formation of a vortex lattice following the insta- 
bility cannot be assumed. Using a simple prediction for 
this frequency, we indicate the regimes in which we ex- 
pect vortex lattice formation to occur. By considering all 
of the key and experimentally tunable quantities in the 
system we outline several accessible routes to generate 
instability and vortex lattices in dipolar condensates. 

This paper is structured as follows. In Section II we 
introduce the mean-field theory and the TF approxima- 
tion for dipolar BECs, in Section III we derive the hy- 
drodynamic equations for a trapped dipolar BEC in the 
rotating frame, and in Section IV we obtain the corre- 
sponding stationary states and discuss their behaviour. 
In Section V we show how to obtain the dynamical sta- 
bility of these states to perturbations, and in Section VI 
we employ the results of the previous sections to discuss 
possible pathways to induce instability in the motion of 
the BEC and discuss the possibility that such instability 
leads to the formation of a vortex lattice. Finally in Sec- 
tion VII we conclude our findings and suggest directions 
for future work. 



II. MEAN-FIELD THEORY OF A DIPOLAR 
BEC 

We consider a BEC with long-range dipolar atomic 
interactions, with the dipolcs aligned in the z direction by 
an external field. The condensate wave function (mean- 
field order parameter) for the condensate tf> = ip(r,t) 
satisfies the GPE which is given by [H |2"51 |2"5] . 



if) 



dif; 



2m 



V 2 + V(r 1 t)+<t> dd (r,t)+gW 



1>, 



(1) 

where m is the atomic mass. The V 2 -term arises from 
kinetic energy and V(r,t) is the external confining po- 
tential. BECs typically feature s-wave atomic interac- 
tions which gives rise to a local cubic nonlinearity with 



coefficient g = AirTi 2 a/m, where a is the s-wave scatter- 
ing length. Note that a, and therefore g, can be experi- 
mentally tuned between positive values (repulsive inter- 
actions) and negative values (attractive interactions) by 
means of a Feshbach resonance [3J [3] . The dipolar inter- 
actions lead to a non-local mean-field potential Q dd (r, t) 
which is given by jS], 



$dd(r, t) = / d 3 r U dd (r - r') p (r' , t) . 



(2) 



where p(r,t) — \ip(r,t)\ 2 is the condensate density and 

l-3cos 2 



U dd (r) 



-in 



(3) 



is the interaction potential of two dipoles separated by a 
vector r, where is the angle between r and the polar- 
ization direction, which we take to be the z-axis. The 
dipolar BECs made to date have featured permanent 
magnetic dipoles. Then, assuming the dipoles to have 
moment d m and be aligned in an external magnetic field 
B — kB, the dipolar coupling is C dd — Hod^ [2H], where 
Po is the permeability of free space. Alternatively, for 
dipolcs induced by a static electric field E = kE, the 
coupling constant C dd — E 2 a 2 /eo [29] [30], where a is 
the static polarizability and eo is the permittivity of free 
space. In both cases, the sign and magnitude of C dd 
can be tuned through the application of a fast-rotating 
external field [31] . 

We will specify the interaction strengths through the 
parameter 



C dd 

V 



(4) 



which is the ratio of the dipolar interactions to the s- 
wave interactions [31] , We take the s-wave interactions 
to be repulsive, g > 0, and so where we discuss negative 
values of e dd , this corresponds to C dd < 0. We will also 
limit our analysis to the regime of —0.5 < e dd < 1, where 
the Thomas-Fermi approach predicts that non-rotating 
stationary solutions are robustly stable [TT]. Outside of 
this regime the situation becomes more complicated since 
the non-rotating system becomes prone to collapse |32j . 

We are concerned with a BEC confined by an elliptical 
harmonic trapping potential of the form, 



V(r) 



1 



[(l-e)x 2 + (l + e)y 2 + j 2 z 2 



(5) 



In the x — y plane the trap has mean trap frequency 
lj±_ and ellipticity e. The trap strength in the axial di- 
rection, and indeed the geometry of the trap itself, is 
specified by the trap ratio 7 = uj z /uj±. When 7 ^> 1 the 
BEC shape will typically be oblate (flattened) while for 
7 <§C 1 it will typically be prolate (elongated), although 
for strong enough dipolar interactions the electrostric- 
tive/magnetostrictive effect can cause a BEC in an oblate 
trap to become prolate itself. 
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The time-dependent GPE can be reduced to 
its time-independent form by making the substitution 
V'fc*) = \J p(n) exp(ifj,t/K), where p is the chemical po- 
tential of the system. We employ the TF approximation 
whereby the kinetic energy of static solutions is taken 
to be negligible in comparison to the potential and in- 
teraction energies. The validity of this approximation in 
dipolar BECs has been discussed elsewhere jTD] . Then, 
the time-independent GPE reduces to, 



V(r) + <S> dd (r) + gp(r) = (J,. 



(6) 



For ease of calculation the dipolar potential Q dd (T.) can 
be expressed as, 



if 



1 



t>dd(r) = -3ge dd I ^(r) + ^p(r) 1 , (7) 
where cf>(r) is a fictitious 'electrostatic' potential defined 

by mm, 



(r) 



1 

47T 



d 3 r'p(r') 



(8) 



This effectively reduces the problem of calculating the 
dipolar potential |2| to the calculation of an electrostatic 
potential of the form Q , for which a much larger theo- 
retical body of literature exists. Exact solutions of Eq. 
(6) for p(r), 0(r) and hence & dd (r) can be obtained for 
any general parabolic trap, as proven in Appendix A of 
Rcf. [12 . In particular, the solutions of p(r) take the 
form 



p(r) 



Po I 1 - ^ - ^ - ) forp(r)>0 (9) 



m 



R 2 



momentum operator. Using this result with the Hamil- 
tonian Hq from Eq. (JTJ we obtain [33J [34] , 



iTi 



dt 



2m 



V 2 + V{r) + $ dd {r,t) + g\*(r,t)f 



0- (x--y- 

i \ dy dx 



*(r,t). 



(11) 



Note that all space coordinates r are those of the rotat- 
ing frame and the time independent trapping potential 
V(r), given by Eq. is stationary in this frame. Mo- 
mentum coordinates, however, are expressed in the lab- 
oratory frame [33J 1311 [35] . 

We can express the condensate mean field in terms 
of a density p(r, i) and phase S(r,t) as ij}(r_,t) — 
\J t) exp[iS(r, t)], and so that the condensate veloc- 
ity is v = (h/m)V_S. Substuting into the time-dependent 



GPE (111 and equating imaginary and real terms leads 



to the following equations of motion, 



dp 

dt 

dv 



= -V • [p (v - x r)] 



-V I -mv ■ v + V(r) + <S> dd {r) 
+ gp — mv ■ \Q X r]) . 



(12) 



(13) 



In the absence of dipolar interactions (& dd = 0) Eqs. ( 12 1 



and ( 13 1 are commonly known as the superfluid hydrody- 
namic equations [6, 7, 8 since they resemble the equation 
of continuity and the Euler equation of motion from dis- 
sipationlcss fluid dynamics. Here we have extended them 
to include dipolar interactions. 

Note that the form of condensate velocity leads to the 
relation, 



where po = 15 N/(8wR x RyR z ) is the central density. Re- 
markably, this is the general inverted parabola density 
profile familiar from the TF limit of non-dipolar BECs. 
An important distinction, however, is that for the dipo- 
lar BEC the aspect ratio of the parabolic solution differs 
from the trap aspect ratio. 



V x v = —V xVS=0, 
m 



(14) 



which immediately reveals that the condensate is irrota- 
tional. The exceptional case is when the velocity poten- 
tial (h/m)S is singular, which arises when a quantized 
vortex occurs in the system. 



III. HYDRODYNAMIC EQUATIONS IN THE 
ROTATING FRAME 



IV. STATIONARY SOLUTION OF THE 
HYDRODYNAMIC EQUATIONS 



Having introduced the TF model of a dipolar BEC we 
now extend this to include rotation and derive hydrody- 
namic equations for the rotating system. We consider 
the rotation to act about the z-axis, described by the ro- 
tation vector £2 where £1 — \Q\ is the rotation frequency 
and the Hamiltonian in the rotating frame is given by, 



H c ft — H — O • L, 



(10) 



where H is the Hamiltonian in absence of the rotation 
and L = —ih(r x V) is the quantum mechanical angular 



We now search for stationary solutions of the hydro- 
dynamic Eqs. (12 1 and (13 1. These states satisfy the 



equilibrium conditions, 



dp 
Of 



= o, 



dv 
di 



= o. 



(15) 



Following the approach of Recati et 
the velocity field ansatz, 

v = aVJxy). 



a. p 



we assume 



(16) 
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Here a is a velocity field amplitude that will provide us 
with a key parameter to parameterise our rotating solu- 
tions. Note that this is the velocity field in the laboratory 
frame expressed in the coordinates of the rotating frame, 
and also that it satisfies the irrotationality condition ( 14 1. 



Combining Eqs. (13 1 and (16 1 we obtain the relation 



M = y (u 2 x x 2 + tfy 2 + w\z 2 ) + gp{r) + <S> dd (r), (17) 
where the effective trap frequencies uj x and Q y are given 

by, 



Cj 2 x = lu 2 j_(1 - e) + a 2 - 2an 



Q 2 V 



wl(l + e) + a 2 + 2afl. 



(18) 



(19) 



The dipolar potential inside an inverted parabola density 
profile ([£]) has been found in Refs. [TSJ G6j to be, 



</r./ 



P0K x Ky 
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x 2 Pwi + y 2 /3 011 + 3z 2 /3 002 
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(20) 



where we have defined the condensate aspect ratios k x = 
R x /R z and n y — R y /R z , and where the coefficients fiijk 
are given by, 



ftijk — 



ds 



[k% + S y +L * ( K 2 + S ) J+S (1 + s ) fe 



i+i 



(21) 



where i, j and fc are integers. Note that for the cylindri- 
cally symmetric case, where k x = n y = k, the integrals 
0ijk evaluate to [55] . 



ftijk 



1F1 (fc+ |,l;i + j + fc+ |;1 - k : 
(l + 2i + 2j + 2k) k 2 ^') 



(22) 



where 2-F1 denotes the Gauss hypergeometric function 
|37j . Thus we can rearrange Eq. (17 1 to obtain an ex- 
pression for the density profile, 



P = 



fj, - f (u 2 x x 2 + Q^y 2 + uj 2 z z 2 ) 

9(1- £dd) 
3ge dd ^g^ [x 2 /3 1Q1 + y 2 f3 011 + 3z 2 /3 002 



RlPooi] 



'(1 -£dd) 



(23) 



Comparing the x 2 , y 2 and z 2 terms in Eq. ([9| and 
Eq. ( 23 1 we find three self-consistency relations that de- 



fine the size and shape of the condensate: 

' lu z \ 2 1 + s d d (§«attj,#L(n - l) 



ft - — - 



c 

u z \ 2 1 + e d d (|«|«iB/3oii - 1) 

c 



2.9PO , 

9 ^ ' 

TTl/, l« 



(24) 

(25) 
(26) 



where £ = 1 — e d d 1 /?oo2 • Furthemore, by in- 
serting Eq. (23) into Eq. (12) we find that stationary 



solutions satisfy the condition, 



(a + fl) (^-^dd ^^ A 



+ (a " fi) ( w= - -£dd 



P01 



c 



(27) 



We can now solve Eq. ( 27 1 to give the velocity field am- 
plitude a for a given e d d, fi and trap geometry. In the 
limit £dd = this amplitude is independent of the s-wave 
interaction strength g and the trap ratio 7. However, in 
the presence of dipolar interactions the velocity field am- 
plitude becomes dependent on both g and 7. For fixed 
e dd and trap geometry, Eq. (27 1 leads to branches of a 
as a function of rotation frequency J7. These branches 
are significantly different between traps that are circular 
(e = 0) or elliptical (e > 0) in the x — y plane, and so 
we will consider each case in turn. Note that we restrict 
our analysis to the range f2 < lu± : for 17 ~ uj± the static 
solutions can disappear, with the condensate becoming 
unstable to a centre-of-mass instability [T5] . 



A. Circular trapping in the x — y plane: e = 

We first consider the case of a trap with no cllipticity 
in the x — y plane (e = 0). In Fig. [T|a) we plot the 
solutions of Eq. ( 27 1 as a function of rotation frequency 
CI for a spherically-symmetric trap 7=1 and for various 
values of e dd . Before dicussing the specific cases, let us 



first point out that for each e dd the solutions have the 
same qualitative structure. Up to some critical rotation 
frequency only one solution exists corresponding to a = 
0. At this critical point the solution bifurcates, giving two 
additional solutions for a > and a < on top of the 
original a = solution. We term this critical frequency 
the bifurcation frequency f2&- 

For e dd — we regain the results of Refs. [15j [16] with 
a bifurcation point at fib = uij_/V2 and, for Q > fit,, 
non-zero solutions given by a — ±y/2fl 2 — uj^jio^ |15j . 
The physical significance of the bifurcation frequency has 
been established for the non-dipolar case and is related to 
the fact that the system becomes energetically unstable 
to the spontaneous excitation of quadrupole modes for 

> Wx/\/2- In the TF limit, a general surface excitation 
with angular momentum Til = TiqiR, where R is the TF 
radius and qi is the quantized wave number, obeys the 
classical dispersion relation tof = (qi / 'm)V rV involving 
the local harmonic potential V = muj\R 2 /2 evaluated at 
R [7] . Consequently, for the non-rotating and non-dipolar 
BEC loi — y/luj±. Meanwhile, inclusion of the rotational 



term in the Hamiltonian (10 1 shifts the mode frequency 
by —Id. Then, in the rotating frame, the frequency of 
the / = 2 quadrupole surface excitation becomes 0)2(^1) = 
y/2u}± — 2f2 [7] . The bifurcation frequency thus coincides 
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FIG. 1: (a) Irrotational velocity field amplitude a of the static 
condensate solutions as a function of the trap rotation fre- 
quency SI in a spherically-symmetric trap (7 = 1 and e = 0). 
Various values of Edd are presented: Edd = —0.49, 0, 0.5 and 
0.99. Insets illustrate the geometry of the condensate in the 
x — y plane, (b) The bifurcation frequency fib (the point at 
which the solutions of a in (a) bifurcate) according to Eq. ( 30 I 
versus trap ratio 7. Plotted are the results for Edd = —0.49, 
-0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 0.9 and 0.99. In (a) and (b) 
Edd increases in the direction of the arrow. 

with the vanishing of the energy of the quadrupolar mode 
in the rotating frame, and the two additional solutions 
arise from excitation of the quadrupole mode for fl > 
wj./V2. 

For the non-dipolar BEC it is noteworthy that fib does 
not depend on the interactions. This feature arises be- 
cause the mode frequencies loi themselves are indepen- 
dent of g. However, in the case of long-range dipo- 
lar interactions the potential <&dd of Eq. (jj gives non- 
local contributions, breaking the simple dependence of 
the force — VV upon R [TT]. Thus we expect the res- 
onant condition for exciting the quadrupolar mode, i.e. 
fib = ui/l (with I = 2), to change with Edd- In Fig. 
[lja) we see that this is the case: as dipole interactions 
are introduced, our solutions change and the bifurcation 
point, fib, moves to lower (higher) frequencies for Edd > 
(sdd < 0). Note that the parabolic solution still satisfies 
the hydrodynamic equations providing —0.5 < Edd < 1- 
Outside of this range the parabolic solution may still ex- 



ist but it is no longer guaranteed to be stable against 
perturbations. 

Density profiles for a = have zero ellipticity in the 
x — y plane. By contrast, the \a\ > solutions have 
an elliptical density profile, even though the trap itself 
has zero ellipticity. This remarkable feature arises due to 
a spontaneous breaking of the axial rotational symme- 
try at the bifurcation point. For a > the condensate is 
elongated in x while for a < it is elongated in y, as illus- 
trated in the insets in Fig. [lja) . In the absence of dipolar 
interactions the |a| > solutions can be intepreted solely 
in terms of the effective trapping frequencies uj x and u> v 
given by Eqs. (fl8| and (Tl9|). The introduction of dipolar 



interactions considerably complicates this picture, since 
they also modify the shape of the solutions. Notably, for 
£dd > the dipolar interactions make the BEC more pro- 
late, i.e., reduce k x and n y , while for Edd < they make 
the BEC more oblate, i.e., increase k x and n y . 

In Fig. [lja) we see that as the dipole interactions are 
increased the bifurcation point fib moves to lower fre- 
quencies. The bifurcation point can be calculated ana- 
lytically as follows. First, we note that for a = the 
condensate is cylindrically symmetric and 
In this case the aspect ratio k is determined by the tran- 
scendental equation [TTJ [T3 121] 



1 



/(«) 
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(£ rfd -l)(« 2 -7 2 ) 

3n 2 £ d d 







where 



/(«) = 



[4 - 3/3, 



000J 



2(1-K 2 ) 



(28) 



(29) 



with /3 oo = (l/v^T?) ln[(l + Vl - « 2 )/(l - x/T^T?)] 
for the prolate case (k < 1), and /?ooo = 
{2/y/ k 2 — 1) arctan[V k 2 — 1] for the oblate case (n > 1). 
For small a — > + , we can calculate the first order correc- 
tions to n x and K y with respect to k from Eqs. ( 24|25 1. 



We can then insert these values in Eq. ( 27 1 and solve for 
$1, noting that in the limit a — ► we have fl — > fib- Thus, 
we find 



n 2 P: 



201 



n b 1 3 2 

w_l V 2 4 l-e dd {l 



101 



|k 2 /3oo2) 



(30) 



In Fig. [ljb) we plot fl b [Eq. ( |3"o| )] as a function of 7 
for various values of Sdd- For Edd = we find that the 
bifurcation point remains unaltered at fib = u! x /y2 as 
7 = uj z /lu x is changed [13 US]- As Edd is increased the 
value of 7 for which fib is a minimum changes from a trap 
shape which is oblate (7 > 1) to prolate (7 < 1). Note 
that for Edd — 0.99 the minimum bifurcation frequency 
occurs at fib ~ 0.55, which is over a 20% deviation from 
the non-dipolar value. For more extreme values of Edd 
we can expect fib to deviate even further, although the 
validity of the inverted parabola TF solution does not 
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FIG. 2: (a) Irrotational velocity field amplitude, a, as a 
function of the trap rotation frequency, £1, for a trap ratio 
7=1 and ellipticity e = 0.025. Various values of Edd are 
presented, Edd = —0.49, 0, 0.5 and 0.99, with E dd increasing 
in the direction of the arrow. Insets illustrate the geometry of 
the condensate in the x — y plane, (b) Backbending point fib 
versus Edd for e = 0.025 and 7 = 0.5 (solid curve), 1.0 (long 
dashed curve) and 2.0 (short dashed curve). 

necessary hold. For a fixed 7 we also find that as Edd 
increases the bifurcation frequency decreases monotoni- 
cally. 



B. Elliptical trapping in the x — y plane: e > 

Consider now the effect of finite ellipticity in the x — y 
plane (e > 0). Rotating elliptical traps have been created 
experimentally with laser and magnetic fields O EJ- 
Following the experiment of Madison et al. [13] we will 
employ a weak trap ellipticity of e = 0.025. In Fig. ^ a) 



we have plotted the solutions to Eq. ( 27 1 for various val- 
ues of Sdd in a 7 = 1 trap. As predicted for non-dipolar 
interactions |15[ 1X6] the solutions become heavily modi- 
fied for e > 0. There exists an upper branch of a > 
solutions which exists over the whole range of fi, and a 
lower branch of a < solutions which back-bends and is 
double- valued. We term the frequency at which the lower 
branch back-bends to be the back-bending frequency fit,. 
The bifurcation frequency in non-elliptical traps can be 



regarded as the limiting case of the back-bending fre- 
quency, with the differing nonclamenture employed to 
emphasise the different structure of the solutions at this 
point. However, for convenience we will employ the same 
parameter for both, fif,. No a — solution exists (for any 
non-zero f2). In the absence of dipolar interactions the 
effect of increasing the trap ellipticity is to increase the 
back-bending frequency fib- Turning on the dipolar in- 
teractions, as in the case of e = 0, reduces Clb for Edd > 0, 
and increases fib for Edd < 0. This is more clearly seen 
in Fig. [2^b) where fl b is plotted versus Edd for various 
values of the trap ratio 7. Also, as in the e = case, 
increasing Edd decreases both k x and fty, i.e. the BEC 
becomes more prolate. 

Importantly, the back-bending of the lower branch can 
introduce an instability. Consider the BEC to be on the 
lower branch at some fixed rotation frequency fl. Now 
consider decreasing Edd- The back-bending frequency fib 
increases and at some point can exceed fl. In other words, 
the static solution of the BEC suddenly disappears and 
the BEC finds itself in an unstable state. We will see 
in Section VI that this type of instability can also be 
induced by variations in 7 and e. 



V. DYNAMICAL STABILITY OF STATIONARY 
SOLUTIONS 



Although the solutions derived above are static solu- 
tions in the rotating frame they are not necessarily stable, 
and so in this section we analyze their dynamical stabil- 
ity. Consider small perturbations in the BEC density and 
phase of the form p = po + 5p and S = So + 8S. Then, 
by linearizing the hydrodynamic equations Eqs. (12 [13 1, 
the dynamics of such perturbations can be described as, 



d 


' 5S~ 




at 


>_ 




where v c 


= v — 



v c ■ V g (1 + EddK) /to 
V-poV [(V • v) + v c ■ V] 



6S 
6p 



(31) 



0. x r and the integral operator K is 



defined as 



(KSp)(r) 



d 2 



Sp(r') dr' 



dz 2 / 4-7T \r ■ 



Sp(r). (32) 



The integral in the above expression is carried out over 
the domain where p Q > 0, that is, the general ellipsoidal 
domain with radii R x , R y , R z of the unperturbed con- 
densate. Extending the integration domain to the region 
where po + Sp > adds higher order effects since it is ex- 
actly in this domain that p — 0(Sp). To investigate the 
stability of the BEC we look for eigenfunctions and eigen- 
values of the operator (31 1: dynamical instability arises 



when one or more eigenvalues A possess a positive real 
part. The size of the real eigenvalues dictates the rate 
at which the instability grows. Note that the imaginary 
eigenvalues of Eq. (31 1 relate to stable collective modes 



of the system [38], e.g. sloshing and breathing, and have 
been analysed elsewhere for dipolar BECs [3S]- In or- 
der to find such eigenfunctions we follow Refs. [03 
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FIG. 3: The maximum positive real eigenvalues of Eq. (311 
(solid curves) for the upper-branch solutions of a as a function 
of fl. We assume e = 0.025, 7 = 1 and N = 3, and present 
various dipolar strengths edd = —0.49, 0, 0.5 and 0.99, with 
Edd increasing in the direction of the arrow. The inset shows 
the full region of dynamical instability in the e — Q, plane for 
£dd = 0. The narrow regions have negligible effect and so 
we only consider the main instability region (bounded by the 
dashed line). 

and consider a polynomial ansatz for the perturbations in 
the coordinates x,y and z, of total degree N. All opera- 
tors in (31 1, acting on polynomials of degree N, result in 



polynomials of (at most) the same degree, including the 
operator K. This latter fact was known to 19th century 
astrophysicists who calculated the gravitational poten- 
tial of a heterogeneous ellipsoid with polynomial density 
[40], |4T] . The integral appearing in Eq. (32 1 is exactly 
equivalent to such a potential. A more recent paper by 
Levin and Muratov summarises these results and presents 
a more manageable expression for the resulting potential 
.42] . Hence, using these results the operator K can be 
evaluated for a general polynomial density perturbation 
bp = x p y q z r , with p, q and r being non-negative integers 
and p + q + r < N. Therefore, the perturbation evolution 
operator (31 ) can be rewritten as a scalar matrix opera- 



tor, acting on vectors of polynomial coefficients, for which 
finding eigenvectors and eigenvalues is a trivial compu- 
tational task. 

Using the above approach we determine the real pos- 
itive eigenvalues of Eq. (31 1 and thereby predict the re- 



gions of dynamical instability of the static solutions. We 
focus on the case of an elliptical trap since this is the ex- 
perimentally relevant case. Recall the general form of the 
branch diagram for this case, i.e. Fig.^a). In the a < 
half-plane, the static solutions nearest the a = axis 
never become dynamically unstable, except for a small 
region ~ due to a centre-of-mass instability of 
the condensate [15]. The other lower branch solutions 
are always dynamically unstable and therefore expected 
to be irrelevant to experiment. Thus, we only consider 
dynamical instability for the upper branch solutions, i.e. 
the branch in the upper half plane (where a > 0). In 
Fig. [3] we plot the maximum positive real eigenvalues of 
the upper branch solutions as a function of £1 for a fixed 
ellipticity e = 0.025. The maximum polynomial pertur- 



bation was set at N = 3, since for this ellipticity it was 
found that higher order perturbations did not alter the 
region of instability, and such modes are therefore not 
displayed. 

For a given Edd and 7 there exists a dynamically un- 
stable region in the e — f2 plane. An illustrative example 
is shown in Fig. pTinset) for Edd = an d 7=1. The in- 
stability region (shaded) consists of a series of crescents 
[16] . Each crescent corresponds to a single value of the 
polynomial degree N, where higher values of N add extra 
crescents from above. At the high frequency end these 
crescents merge to form a main region of instability, char- 
acterised by large eigenvalues. At the low frequency end 
the crescents become vanishingly thin and are charac- 
terised by very small eigenvalues which are at least one 
order of magnitude smaller than in the main instabil- 
ity region [19J . As such these regions will only induce 
instability in the condensate if they are traversed very 
slowly. This was confirmed by numerical simulations in 
Ref. [T3] where it was shown that the narrow instability 
regions have negligible effect when ramping f2 at rates 
greater than dfl/dt — 2 x 10~ 4 or^. It is unlikely that 
an experiment could be sufficiently long-lived for these 
narrow instability regions to play a role. For this reason 
we will subsequently ignore the narrow regions of insta- 
bility and define our instability region to be the main 
region, as bounded by the dashed line in Fig. [3pnset). 
For the experimentally relevant trap ellipticities e^O.l 
the unstable region is defined solely by the N = 3 per- 
turbations. 

We define the lower bound of the instability region to 
be 0^ (this corresponds to the dashed line in the inset). 
This is the key parameter to characterise the dynamical 
instability. As Edd is increased f2j decreases and, accord- 
ingly, the unstable range of widens. Note that the 
upper bound of the instability region is defined by the 
endpoint of the upper branch at ~ luj_ . 



VI. ROUTES TO INSTABILITY AND VORTEX 
LATTICE FORMATION 

A. Procedures to induce instability 

For a non-dipolar BEC the static solutions and their 
stability in the rotating frame depend only on rotation 
frequency fl and trap ellipticity e. Adiabatic changes in e 
and f2 can be employed to evolve the condensate through 
the static solutions and reach a point of instability. In- 
deed, this has been realized both experimentally [T51 IT4] 
and numerically [T7] [H], with excellent agreement to the 
hydrodynamic predictions. For the case of a dipolar BEC 
we have shown in Sections IIVI and [V] that the static so- 
lutions and their instability depend additionally on the 
trap ratio 7 and the interaction parameter Edd- Since 
all of these parameters can be experimentally tuned in 
time, one can realistically consider each parameter as a 
distinct route to traverse the parameter space of solutions 
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and induce instability in the system. 

Examples of these routes are presented in Fig. [4] 
Specifically, Fig. [3] shows the static solutions a of Eq. 
J27} as a function of [Fig. ga)] , e [Fig. gfb)], e dd [Fig. 
gfcj] and 7 [Fig.gd)]. In each case the remaining three 
parameters are fixed at e = 0.025, 7 = 1, f2 = 0.7w_i_, 
and £dd — 0.99. Dynamically unstable solutions are indi- 
cated with red circles. Grey arrows mark routes towards 
instability (the point of onset of instability being marked 
by an asterisk), where the free parameter f2, e, Edd, or 7 
is varied adiabatically until cither a dynamical instabil- 
ity is reached, or the solution branch backbends and so 
ceases to exist. For solutions with a > 0, the instability 
is always due to the system becoming dynamically unsta- 
ble (dashed arrows), whereas for a < the instability is 
always due to the solution branch backbending on itself 
(solid arrows) and so ceasing to exist. Numerical studies 
[T8"] indicate that these two types of instability involve 
different dynamics and possibly have distinct experimen- 
tal signatures. 

Below we describe adiabatic variation of each parame- 
ter in more general detail, beginning with the established 
routes towards instability in which (i) f2 and (ii) e are 
varied, and then novel routes based on adiabatic changes 
in (iii) Edd and (iv) 7. In each case it is crucial to con- 
sider the behaviour of the points of instability, namely 
the back-bending point fib and the onset of dynamical 
instability of the upper branch 

i) Adiabatic introduction of f2: The relevant pa- 
rameter space of e and f2 is presented in Fig. [5ja,) , with 
the instability frequencies fif,(e) (solid curves) and fii(e) 
(dashed curves) indicated. For a BEC initially confined 
to a non-rotating trap with finite ellipticity e, as the ro- 
tation frequency f2 is increased adiabatically the BEC 
follows the upper branch solution [Fig. ga) (dashed ar- 
row)]. This particular route traces out a horizontal path 
in Fig. [5 a) until it reaches ^i(e), where the stationary 
solution becomes dynamically unstable. For the specific 
parameters of Fig. ga) the system becomes unstable at 

= f2j(e) w 0.65cjj^. More generally, Fig. ga) shows 
that as Edd is increased, f2j(e) is decreased and as such 
instabilities in the stationary solutions will occur at lower 
rotation frequencies. At e ~ 0.1, the curve for displays 
a sharp kink, arising from the shape of the dynamically 
unstable region, as shown in Fig. g inset). 

ii) Adiabatic introduction of e: Here we begin with 
a cylindrically symmetric (e = 0) trap, rotating at a fixed 
frequency f2. The trap ellipticity e is then increased adi- 
abatically and in the phase diagram of Fig.ga) the BEC 
traces out a vertical path starting at e — 0. The ensuing 
dynamics depend on the trap rotation speed relative to 
n b (e = 0): 

(a) For < f2b(e = 0) the condensate follows the upper 
branch of the static solutions shown in Fig. ga). This 
branch moves progressively to larger a. For < fij(e) 
the BEC remains stable but as e is increased further 
the condensate eventually becomes dynamically unsta- 
ble. Figure g a) shows that as Edd is increased n$(e) is 
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FIG. 4: Stationary states in the rotating trap characterised 
by the velocity field amplitude a, determined from Eq. (27 1. 



Dynamically unstable solutions are marked with red circles. 
In each of the figures the trap rotation frequency Q (a), trap 
ellipticity e (b) , dipolar interaction strength Edd (c) and axial 
trapping strength 7 (d) are varied adiabatically, whilst the 
remaining parameters remain fixed at Q — 0.7luj_, e = 0.025, 
£dd = 0.99, and 7=1. The adiabatic pathways to instability 
(onset marked by red asterisk) are schematically shown by 
the dashed and solid arrows. Dashed arrows indicate a route 
towards dynamical instability, whereas solid arrows indicate 
an instability due to disappearance of the stationary state. 



decreased and as such the dynamical instability of the 
stationary solutions occurs at a lower trap ellipticity. 

(b) For > flb(e = 0) the condensate accesses the 
lower branch solutions nearest the a = axis. These so- 
lutions are always dynamically stable and the criteria for 
instability is instead determined by whether the solution 
exists. As e is increased the back-bending frequency fife(e) 
increases. Therefore, when e exceeds some critical value 
the lower branch solutions disappear for the chosen value 
of rotation frequency f2. This occurs when < f2f,(e). 
Figure g a) shows that as Edd is increased fi(,(e) is de- 
creased and as such instabilities in the system will occur 
at a higher trap ellipticity. At this point the parabolic 
condensate density profile no longer represents a stable 
solution. The particular route indicated in Fig. gb) is 
included in Fig. ga) as a vertical, solid grey arrow. 

iii) Adiabatic change of Edd- The relevant param- 
eter space of Edd and Q is shown in Fig. 5(b) for several 
different trap ratios. Consider that we begin from an ini- 
tial BEC in a trap with finite ellipticity e = 0.025 and 
rotation frequency £1. (This can be achieved, for exam- 
ple, by increasing e from zero at fixed f2.) Then, by 
changing £d d adiabatically an instability can be induced 
in two ways: 

(a) For ft < flb(edd) the condensate follows the up- 
per branch solutions until they become unstable. This 
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FIG. 5: (a) Phase diagram of e versus Qb (solid curves) and 
Qi (dashed curves) for 7=1 and (i) Edd = —0.49, (ii) 0.5 and 
(iii) 0.99. (b) Phase diagram of Edd versus Clb (solid curves) 
and fij (dashed curves) for e = 0.025 and (i) 7 = 0.5, (ii) 1 
and (iii) 2. (c) Phase diagram of 7 versus fit (solid curves) 
and fli (dashed curves) for e = 0.025 and (i) Edd = —0.49, 
(ii) 0, (iii) 0.5 and (iv) 0.99. In each case the solid (dashed) 
arrows depict the routes to instability shown in Fig. 4. 

route to instability in shown in Fig. 4(c) by the dashed 
arrow, with the corresponding path in Fig. 5(b) shown 
by the vertical dashed arrow. Thus for CI < Cli(Edd) the 
motion remains stable. However, for CI > Cli(sdd) the up- 
per branch becomes dynamically unstable. In Fig. [5jb) 
Qi(£dd) (dashed curves) is plotted for different trap ra- 



tios. As can be seen, the stable region of the upper branch 
becomes smaller as Edd is increased. 

(b) For CI > Clb(£dd) the condensate follows the lower 
branch solutions nearest the a = axis. These so- 
lutions are always stable and hence an instability can 
only be induced when this solution no longer exists, i.e. 
Cl < Qb( s dd)- Figure [5|b) shows Cl b (Edd) (solid curves) 
for various trap aspect ratios. As can be seen the back- 
bending frequency Clb decreases as Edd is increased. Thus 
if Edd is increased the system will remain stable. However 
if Edd is decreased then the system will become unstable 
when Cl = Cl b (E d d)- 

iv) Adiabatic change of 7: Figure [5jc) shows the 
parameter space of 7 and Cl. Consider, again, an ini- 
tial stable condensate with finite trap rotation frequency 
Cl and ellipticity e = 0.025. Then through adiabatic 
changes in 7 the condensate can traverse the parame- 
ter space and, depending on the initial conditions, the 
instability can arise in two ways: 

(a) For Cl < flb(j) the condensate exists on the up- 
per branch. It is then relevant to consider the onset of 
dynamical instability ^(7) (dashed curves in Fig. |5jc)). 
Providing Cl < ^(7) the solution remains dynamically 
stable. However, once Cl > ^(7) the upper branch solu- 
tions become unstable. 

(b) For Cl > Clbij) the condensate exists on the lower 
branch nearest the a — axis. These solutions are always 
dynamically stable and instability can only occur when 
the motion of the back-bending point causes the solution 
to disappear. This occurs when CI < ^(7), with ^(7) 
shown in Fig. [5^c) by solid curves for various dipolar 
interaction strengths. 

These two paths to instability are shown in Fig. J4jd) 
and are also indicated in Fig.[5jc) as vertical grey arrows, 
where the dashed (solid) arrow corresponds to the a > 
(a < 0) path. 



B. Is the final state of the system a vortex lattice? 

Having revealed the points at which a rotating dipo- 
lar condensate becomes unstable we will now address the 
question of whether this instability leads to a vortex lat- 
tice. First, let us review the situation for a non-dipolar 
BEC. The presence of vortices in the system becomes en- 
ergetically favorable when the rotation frequency exceeds 
a critical frequency Cl v . Working in the TF limit, with 
the background density taking the parabolic form (|9| , f2„ 
can be approximated as |49j . 



Ch<i; 



2 mR 2 



h , 0.67i? 
In ■ 



(33) 



Here the condensate is assumed to be circularly symmet- 
ric with radius R, and £ s = Ti/^/2mpog is the healing 
length that characterises the size of the vortex core. For 
typical condensate parameters Cl v ~ 0.4cj^. It is ob- 
served experimentally, however, that vortex lattice for- 
mation occurs at considerably higher frequencies, typi- 
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FIG. 6: The relation between the instability frequencies, Q,b 
(long dashed red curve) and fij (short dashed curve) , and the 
critical rotation frequency for vorticity Q v (solid curve) for 
(a) an oblate trap y = 10 and (b) a prolate trap 7 = 0.1. The 
instability frequencies are based on a trap with ellipticity e = 
0.025 while Q v is obtained from Eq. ( 33 1 under the assumption 
of a 52 Cr BEC with 150, 000 atoms and scattering length a s = 
5.1nm in a circularly symmetric trap with lj± = 2n x 200Hz. 



cally ~ 0.7lu±. This difference arises because above Sl v 
the vortex-free solutions remain remarkably stable. It is 
only once a hydrodynamic instability occurs (which oc- 
curs in the locality of f2 « 0.7oj±) that the condensate has 
a mechanism to deviate from the vortex-free solution and 
relax into a vortex lattice. Another way of visualising this 
is as follows. Above fl v the vortex-free condensate resides 
in some local energy minimum, while the global mini- 
mum represents a vortex or vortex lattice state. Since 
the vortex is a topological defect, there typically exists 
a considerable energy barrier for a vortex to enter the 
system. However, the hydrodynamic instabilities offer a 
route to navigate the BEC out of the vortex-free local 
energy minimum towards the vortex lattice state. 

Note that vortex lattice formation occurs via non- 
trivial dynamics. The initial hydrodynamic instability in 
the vortex-free state that we have discussed in this paper 
is only the first step [18 . For example, if the condensate 
is on the upper branch of hydrodynamic solutions (e.g. 
under adiabatic introduction of fi) and undergoes a dy- 
namical instability, this leads to the exponential growth 
of surface ripples in the condensate [T31 [13]. Alterna- 
tively, if the condensate is on the lower branch and the 



static solutions disappear (e.g. following the introduc- 
tion of e) the condensate undergoes large and dramatic 
shape oscillations. In both cases the destabilisation of 
the vortex-free condensate leads to the nucleation of vor- 
tices into the system. A transient turbulent state of vor- 
tices and density perturbations then forms, which subse- 
quently relaxes into a vortex lattice configuration [T5J H5J . 

In the presence of dipolar interactions, however, the 
critical frequency for a vortex depends crucially on the 
trap geometry 7 and the strength of the dipolar interac- 
tions Edd- Following reference |21j we will make a sim- 
ple and approximated extension of Eq. ( 33 1 to a dipo- 
lar BEC. We will consider a circularly-symmetric dipo- 
lar condensate with radius R — R x — R y that satisfies 
Eqs. ( 24 )-( 26 1, and insert this into Eq. (33 1 for the con- 



densate radius. This method still assumes that the size of 
the vortex is characterised by the s-wave healing length 
£ s . Although one does expect the dipolar interactions to 
modify the size of the vortex core, it should be noted that 
Eq. ( 33 1 only has logarithmic accuracy and is relatively 
insensitive to the choice of vortex core length scale. The 



dominant effect of the dipolar interactions in Eq. ( 33 1 
comes from the radial size and is accounted for. Note also 
that this expression is for a circularly symmetric system 
while we are largely concerned with elliptical traps. How- 
ever we will employ a very weak ellipticity e = 0.025 for 
which we expect the correction to the critical frequency 
to be correspondingly small. 

As an example, we take the parameter space of rotation 
frequency il and dipolar interactions edd- We first con- 
sider the behavior in a quite oblate trap with 7=10. In 
Fig.[6ja) we plot the instability frequencies fli and flj for 
this system as a function of the dipolar interactions Edd- 
Depending on the specifics of how this parameter space 
is traversed, either by adiabatic changes in f2 (vertical 
path) or Edd (horizontal path), the condensate will be- 
come unstable when it reaches one of the instability lines 
(short and long dashed lines). These points of instability 
decrease weakly with dipolar interactions and have the 
approximate value fli ss f2f, ss 0.75wj_. On the same plot 
we present the critical rotation frequency Q v according 
to Eq. ( 33 1 . In order to calculate this we have assumed a 
BEC of 150,000 52 Cr atoms confined within a trap with 
w± = 2tt x 200Hz. In this oblate system we see that the 
dipolar interactions lead to a decrease in f2 OJ as noted in 
[2Tj . This dependence is very weak at this value of 7, 
and throughout the range of Sdd presented it maintains 
the approximate value £l v ~ 0.1wj_. Importantly these 
results show that when the condensate becomes unstable 
a vortex/vortex lattice state is energetically favored. As 
such, we expect that in an oblate dipolar BEC a vortex 
lattice will ultimately form when these instabilities are 
reached. 

In Fig. |6jb) , we make a similar plot but for a prolate 
trap with 7 = 0.1. The instability frequencies show a 
somewhat similar behaviour to the oblate case. How- 
ever, Q v is drastically different, increasing significantly 
with Edd- We find that this qualitative behaviour occurs 



11 



consistently in prolate systems, as noted in [21] ■ This 
introduces two regimes depending on the dipolar inter- 
actions. For e dd <0.8, O iib > Q v , and so we expect a 
vortex/vortex lattice state to form following the insta- 
bility. However, for £dd^0.8 we find an intriguing new 
regime in which f^./, < fi„. In other words, while the in- 
stability in the vortex-free parabolic density profile still 
occurs, a vortex state is not energetically favorable. The 
final state of the system is therefore not clear. Given that 
a prolate dipolar BEC is dominated by attractive inter- 
actions (since the dipoles lie predominantly in an attrac- 
tive end-to-end configuration) one might expect similar 
behavior to the case of conventional BECs with attrac- 
tive interactions (g < 0) where the formation of a vor- 
tex lattice can also be energetically unfavorable. Sugges- 
tions for final state of the condensate in this case include 
centre-of-mass motion and collective oscillations, such as 
quadrupole modes or higher angular momentum-carrying 
shape excitations [IS 35J 137] . However the nature of the 
true final state in this case is beyond the scope of this 
work and warrants further investigation. 

VII. CONCLUSIONS 

By calculating the static hydrodynamic solutions of 
a rotating dipolar BEC and studying their stability, we 
have predicted the regimes of stable and unstable motion. 
In general we find that the backbcnding or bifurcation 
frequency decreases with increasing dipolar interac- 
tions. In addition, the onset of dynamical instability in 
the upper branch solutions, f2j, decreases with increas- 



ing dipolar interactions. Furthermore these frequencies 
depend on the aspect ratio of the trap. 

By utilising the novel features of dipolar condensates 
we detail several routes to traverse the parameter space 
of static solutions and reach a point of instability. This 
can be achieved through adiabatic changes in trap rota- 
tion frequency £1, trap ellipticity e, dipolar interactions 
Edd and trap aspect ratio 7, all of which are experimen- 
tally tunable quantities. While the former two methods 
have been employed for non-dipolar BECs, the latter two 
methods are unique to dipolar BECs. In an experiment 
the latter instabilities would therefore demonstrate the 
special role played by dipolar interactions. Furthermore, 
unlike for conventional BECs with repulsive interactions, 
the formation of a vortex lattice following a hydrody- 
namic instability is not always favored and depends sen- 
sitively on the shape of the system. For a prolate BEC 
with strong dipolar interactions there exists a regime in 
which the rotating spheroidal parabolic Thomas-Fermi 
density profile is unstable and yet it is energetically un- 
favorable to form a lattice. Other outcomes may then 
develop, such as a centre-of-mass motion of the system 
or collective modes with angular momentum. However, 
for oblate dipolar condensates, as well as prolate con- 
densates with weak dipolar interactions, the presence of 
vortices is energetically favored at the point of instability 
and we expect the instability to lead to the formation of 
a vortex lattice. 
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(NGP) and the Natural Sciences and Engineering Re- 
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